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Abstract 



We study the dynamics of relaxation and thermalization in an exactly 
solvable model of a particle interacting with an harmonic oscillator bath. 
Our goal is to understand the effects of non-Markovian processes on the re- 
laxational dynamics and to compare the exact evolution of the distribution 
function with approximate Markovian and Non-Markovian quantum kinetics. 
There are two different cases that are studied in detail: i) a quasiparticle (reso- 
nance) when the renormalized frequency of the particle is above the frequency 
threshold of the bath and ii) a stable renormalized 'particle' state below this 
threshold. The time evolution of the occupation number for the particle is 
evaluated exactly using different approaches that yield to complementary in- 
sights. The exact solution allows us to investigate the concept of the formation 
time of a quasiparticle and to study the difference between the relaxation of 
the distribution of bare particles and that of quasiparticles. For the case of 
quasiparticles, the exact occupation number asymptotically tends to a statisti- 
cal equilibrium distribution that differs from a simple Bose-Einstein form as a 
result of off-shell processes whereas in the stable particle case, the distribution 
of particles does not thermalize with the bath. We derive a non-Markovian 
quantum kinetic equation which resums the perturbative series and includes 
off-shell effects. A Markovian approximation that includes off-shell contri- 
butions and the usual Boltzmann equation (energy conserving) are obtained 
from the quantum kinetic equation in the limit of wide separation of time 
scales upon different coarse-graining assumptions. The relaxational dynamics 
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predicted by the non-Markovian, Markovian and Boltzmann approximations 
are compared to the exact result. The Boltzmann approach is seen to fail in 
the case of wide resonances and when threshold and renormalization effects 
are important. 

76.20.+q;72.10.Bg;72.15.Lh 



2 



Typeset using REVTgX 



I. INTRODUCTION AND MOTIVATION 



Recent advances in semiconductor femtosecond spectroscopy [|I]J2j highlight the need 
for a deeper theoretical understanding of the relaxational dynamics of hot carriers that 
goes beyond Boltzmann kinetics. Boltzmann or semiconductor Bloch equations, are based 
on strict energy conservation and result in a Markovian description as a consequence of 
averaging over microscopic time scales. On short time scales, the time-energy uncertainty 
principle comes into play and off-shell (non-energy conserving) processes lead to quantum 
kinetic equations with memory effects, i.e. non-Markovian effects. 

Ultrafast relaxation in semiconductors are typically studied by exciting a semiconductor 
sample with a femtosecond laser ||- [|J. The subsequent dynamics of the photoexcited 
carriers is then studied by measuring the optical or transport properties of the sample at 
different time delays. These experiments demonstrate the breakdown of Boltzmann kinetics 
for periods less than the optical lattice oscillation period (around 115 fs in GaAs 0) and 
emphasize the need for a quantum kinetic description of the relaxational dynamics. 

Motivated by these new developments, there is a rekindled interest on a deeper theoretical 
understanding of quantum kinetics and critical analysis of transport and kinetic approaches 
are beginning to emerge f~ 
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In particular recently the initial stages of pre-equilibration 
during which quasiparticle correlations begin to build had been investigated in a many body 
system ||. The pre-equilibrium stage cannot be studied within a Boltzmann approach 
because the early time dynamics depends on the initial preparation of the state and is 
determined by virtual processes that do not conserve energy on short times (off-shell). 

Besides semiconductor systems, the interest in quantum kinetics is truly interdisciplinary: 
in dense plasma ||, nuclear matter []8,10,11] and high energy physics and cosmology fl2| - 
| 14| to cite but a few applications. 

A very powerful method to derive quantum kinetic equations uses the non-equilibrium 
Green's functions within the Keldysh formalism |L5| which leads to the Kadanoff-Baym 
equations fllE| . In order to derive quantum kinetic equations some assumptions must be 
invoked. Usually the generalized Kadanoff-Baym ansatz |17|- [21] with renormalized one 
particle Green's function propagators |£| is used to relate the two-time correlation functions 
with the one-time distribution function. An alternative approach to derive quantum kinetic 
equations is by truncating the Bogolyubov-Born-Green-Kirkwood-Yvon (BBGKY) hierarchy 



Although there is ample experimental and numerical confirmation of kinetics described 
by Boltzmann or kinetic Bloch equations in processes in which there is a wide separation 
between relaxational and microscopic time scales, the situation for non-Markovian quantum 
kinetics is less well understood. 

An important limitation in the numerical study of non-Markovian kinetic equations is 
the very intensive computational requirements to analyze integro-differential equations with 
Thus it is important to try to test relaxation via non-Markovian quantum ki- 
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memory 

netic equations in systems which afford an exact solution. Recently non-Markovian quantum 
kinetics has been studied within the context of hot electron relaxation ][H|] in one dimension. 
This model affords an exact solution via bosonization and allows a direct comparison to an 
approximate kinetic treatment. Furthermore, improved transport equations that include the 
effects of correlations leading to a non-Markovian description have been recently proposed 
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24J and compared to available exact solutions in low dimensional models. 

Thus the current experimental efforts in femtosecond relaxation in semiconductors and 
the necessity for a deeper understanding of quantum kinetics via non-Markovian transport 
equations justifies the study of model systems that can be solved exactly and thus provide 
a testing ground for the different types of approximations. 

The goal of this article is to study the description of thermalization and relaxational 
dynamics in a simple and exactly solvable many body theory to obtain a deeper under- 
standing of the off-shell processes (not energy conserving on short time scales) involved in 
thermalization and to provide a yardstick to test different approximations. 

The aspects that we seek to study in this article are the following: 

• How do off-shell effects modify the dynamics of thermalization and relaxation? By off- 
shell we here refer to processes that do not conserve energy on short time scales and 
threshold effects that are not incorporated in the usual Boltzmann equation. These 
are responsible for quasiparticle properties such as widths and wave function renor- 
malization. 

• A detailed understanding of the relaxation of quasiparticles vs. that of bare and 
dressed particles and to explore the definition of a quasiparticle distribution function 
that is valid beyond the narrow width approximation. The model under consideration 
also allows us to study the formation time of the quasiparticle. 

• A comparison of the validity of Markovian (coarse grained) approximations including 
the Boltzmann equation, to a non-Markovian description of relaxation which is a 
simplified form of the Kadanoff-Baym equations. 

Although we anticipate that the answer to many of these questions will in general depend 
on the details of the microscopic model, we propose to study a model of a particle (harmonic 
oscillator) interacting linearly with a bath of harmonic oscillators. As it will be seen in what 
follows this model bears many of the properties of more realistic interacting systems. By 
studying different couplings between the particle and the bath, we provide answers to these 
(and other) questions and obtain further intuition into more complex situations. 

In section II we introduce the model and discuss the different approaches to study the 
dynamical evolution of the distribution function. In section III we analyze the dynamics from 
the point of view of the time evolution of an initially prepared density matrix which allows 
us to establish contact with the fluctuation dissipation relation. Here we distinguish between 
bare particle and dressed particle and quasi-particle distributions. The exact solution of the 
Heisenberg equations of motion is presented in section IV. In section V we study in detail 
the long time dynamics of the distribution function. In section VI we discuss the exact 
solution in terms of normal modes and analyze the definition of the quasiparticle number 
operator that describes the relaxational dynamics. In section VII we analyze the approximate 
relaxational dynamics in terms of i) the Boltzmann equation, ii) the non-Markovian quantum 
kinetic equation and iii) a Markovian approximation to the quantum kinetic equation. We 
provide a numerical comparison of the exact and approximate kinetics in section VIII. Our 
conclusions are summarized in section IX. 
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II. THE MODEL 



As stated in the introduction, we seek to study aspects of quantum kinetics in a model 
that allows to compare approximate treatments of the relaxational dynamics to exact solu- 
tions. The model that we choose to describe this situation is that of an oscillator of bare 
frequency uq (representing the physical mass of the in particle states before the interaction) 
coupled linearly with a bath with an infinite number of degrees of freedom given by harmonic 
oscillators with frequencies Uk- Although this is a drastic simplification of microscopic inter- 
acting theories, this model continues to serve as a testing ground for studies of dissipation 



in quantum systems 26 



The Lagrangian is given by 

L[q,Q k ] = \(? - <4f) + \Y,(qI-"1 Qt) Q k , 

k k 

where the different coefficients of q 2 ; Q\ (oscillator masses) had been absorbed by a canonical 
transformation into a redefinition of the couplings C^. We now refer to the oscillator q as 
the 'system' i.e. the degree of freedom whose dynamics we are interested in studying, and 
the oscillators Qk as the 'bath', these will be integrated out in the non-equilibrium effective 
action. This model also describes the interaction of an electron with a phonon or photon 



bath in the dipole approximation f3T| , |32 



We will eventually take the limit in which the bath oscillators are distributed continuously 
by introducing the bath spectral density J(ui) and where appropriate replacing the discrete 
distribution with a continuum one in the following manner: 



in such a way that 



7T C 2 

J M = o £ — - ^) 

1 k Uk 



J2C 2 k f(u k )^- [ du u J(u) f{u>). (2.1) 

k 71 J 



Our main goal is to study the evolution of the number of excitations or 'particle distri- 
bution' associated with the quanta of the system. Anticipating self-energy renormalization 
effects by the bath, we define a reference frequency Q and introduce the operator that counts 
the number of quanta of the system's degrees of freedom associated with this frequency 
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1 



n{t) = ^ P z (t) + OV(f) - - , (2.2) 



2 



where p(t) is the momentum of the particle. The reference frequency could either be taken 
to be the bare frequency Uq, or the frequency renormalized by the interaction with the bath, 
we will leave this choice unspecified for the moment. 

Since the theory is quadratic we can resort to a number of different ways to study the 
dynamical evolution: 

1. Given an initial density matrix we can evolve it in time exactly and obtain all of the 
non-equilibrium correlation functions. 
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2. The Heisenberg equations of motion for the operators can be solved exactly and again 
we can obtain any correlation function. 

3. The normal modes can be found exactly, from which we can find the exact ground 
state and also obtain the operators that create the particle or quasiparticle states to 
study the asymptotic evolution of non-equilibrium states. 

4. We compute exactly the expectation value of the proposed number operator in the 
canonical ensemble of the system plus bath and compare the result to the asymp- 
totic form of the non-equilibrium distribution function. This allows an unequivocal 
description of thermalization in terms of the density matrix. 

We will pursue all of the above different approaches, since each particular method pro- 
vides different insights and the main goal is to understand this simpler model in detail to 
provide intuition into more realistic cases. 



III. TIME EVOLUTION OF AN INITIAL DENSITY MATRIX 



The first method is to calculate the time evolution of the reduced density matrix, p r (t), 
of the particle that has been prepared at some initial time £j. 

This can be achieved by treating the infinite set of harmonic oscillators, Qk, as a 'bath' 
and obtaining an influence functional p7| - pl[] by tracing out the bath degrees of freedom. 
We assume that the total density matrix for the particle-bath system decouples at the initial 



time ti, i.e. 



p(U) = p s (ti) <g> p R (ti), 



where pii(ti) is the density matrix of the bath which describes infinite set of harmonic 
oscillators in thermal equilibrium at a temperature T and p s (ti) is the density matrix of 
the particle which is taken to be that of a harmonic oscillator in thermal equilibrium at 
temperature T . More complicated initial density matrices, including correlations between 



system and bath degrees of freedom can be studied by following the methods found in [28 



The complete information of non-equilibrium expectation values and correlation func- 
tions is completely contained in the time dependent density matrix 

P (t) = u(t,u) P {u) u- l {t,u) 

with U(t,ti) the time evolution operator. Real time non-equilibrium expectation values and 
correlation functions can be obtained via functional derivatives with respect to sources of 
the generating functional [ |i~5| , [33|l 



Tr 



f/(oo,t i; j + ) p{U) U-^oOtUJ-)] /Trp(ti) 



where j are sources coupled to the particle coordinate. This generating functional is readily 
obtained using the Schwinger-Keldysh method which involves a path integral in a complex 
contour in time Hl5|,33|. Real time, non-equilibrium Green's functions are now obtained as 



functional derivatives with respect to the sources. There are four types of free propagators 
MM 



6 



Qt(t)Qt(t')) = -*Gt + (t,t>) = -i[G>(t,t')9(t-t') + G<(t,t')e(t'-t)} 

QkWQkW) = -iGTW) = -i[G^(t,t')e(t'-t) + G<(t,t')9(t-t')} 

QtWQkW) = iG+-(t,t') = -iG<(t,t') 

QkWQtW) = iG^ + (t,t>) = -iG>(t,t') = -iG<(t',t), 



(3.1) 



where the signs ± in the above expressions correspond to the fields and sources on the 
forward (+) and backward (— ) branches and 



G>(t,t') 
G<(t,t') 



2uj k 
% 

2u) k 



[l + N k ) exp {-iuj k (t - t')} + N k exp{iaj k (t-t')} 

[I + N k ) exp {ioo k (t - t')} + N k exp{-iuj k (t-t')} 
1 



exp{(3u k } - 1 



(3.2) 



A. The reduced density matrix 



The reduced density matrix, p r (t), is defined as [pTfl - [31 

Tr R p(t) 



p r (t) 



Trp R (ti 



where the subscript R in Tr R refers to tracing over the bath degrees of freedom. Taking the 
trace over Q k , one obtains the reduced density matrix in terms of the influence functional 
0-0, ^[g+g-] 



•[q,q';t] = J dq x dq 2 po[?i,?2] J Vq + Vq exp li J dt (^L [q + ] - L [q ])| x J"[g + ,g" 



(3.3) 



with the following boundary conditions on the fields: q + (ti) = q\\ q + (oo) — q; q (oo) 
q'\ q~(ti) = q 2 - Po[qi,q2] is the initial density matrix of the particle and 



T[q\q-] = expi^C, 2 / dt J dt> £ q\t)Gf (t, t)q\t>) \ 

y k a,b ) 



; a, b 



1 r 



We will choose the initial density matrix of the particle to be that of an harmonic 
oscillator of reference frequency Q in thermal equilibrium at temperature T given by 



Po [91,92] 



1 



V 2ttct 
x exp 



exp {ipi(qi - q 2 )} 
{ Q 



2sinh [(3 Q,) 



((li - Qi) 2 + (<?2 - Qi) 2 ) cosh [f3 Q] - 2{q l - ^)(g 2 - qi 



where g« and pi are respectively the average position and momentum of the particle, (3q 
1/T and 



a 



1 



coth 



Pott 



1 + 2n(0) 
2SL 



n(0) 



1 



e /3 ^ _ 1 ' 



The reference frequency Q will allow to understand the different features of the dynamics 
of the dressed particle in the medium, rather than the bare particle with frequency uq. We 
will specify this reference frequency below when we study the dynamics in detail. 



(3.4) 



Using the Wigner coordinates ||27|| - pi ] which are defined as 

*(0 = \ (<? + (0 + <r(0) . KO = q + {f) - <f (0, 



the integrals in eq. (3.3) can be evaluted easily and one obtains the reduced density matrix 



Pr[x f ,rf,t\ 



: exp 



-g(t) B(t) 



a , x B 2 (t) 
— + R-(t) U- 



(g-(t)) 



X f Tf + I 



2 Ait) 



~X 4 



B(t) 



git) 2A{t)\~>-> ' ^^Ait)^^^ 
\f^{pi9{t) + Qi9(t)) x f - jjfrj(Pi9(t) + 9i9 



4A(tf f 



2A(t) 



where 



A(t) 
B{t) 



^V(t) + \ R ++it) + °-g\t) 



(3.5) 



R++(t) = C dt! C dfgit - t') Kit' - t") git - t' 
Jo Jo 

R-(t) = f dt' f dt" €^1 Kit' - t") ^ 
Jo Jo g it) g it) 

g-{t") 
g-it) 

COS [Uk(t' — t")] 



R + ~(t)= f dt' [ dt" g(t-t')K(t' -t 
Jo Jo 

(3uj k 



Kit' - 1" 

g-it' 



coth 



? 2^ 

git)gjt-t') - g{t)g{t-f) 
g(t)g(t) ~ 9 2 (t) 



(3.6) 



The dynamics of the reduced density matrix is completely determined by the function g{t) 
which satisfies the following differential equation 



git)+^og(t) ~ f dfT.it -t')g{t') = 
Jo 



(3.7) 



with initial conditions g{0) = g(0) = and g{0) = 1. The kernel £(£ — t') is the retarded 
self energy of the system degree of freedom and it is given by 



S 



s(t - 1) = e(t - E — sin - 0] • ( 3 - 8 ) 

We will postpone the computation of the function g(t) to the next section and we will 
specify the spectral density for the bath in a later section wherein we will compare exact 
results to different approximations. 

Having obtained the reduced density matrix, we can now obtain the expectation values 
of q 2 (t) and p 2 (t) and to compute the expectation value of the number operator (|2.2j) , which 
after some straightforward algebra is shown to be given by 

(„(()> =-i + K(t) + §(rf + S2V) J 2 W + (j Y^!)Ht) 

+ P?+ " 2 2 g +2g) 9 ! « + *f (SW + tf 9(i) ) (3.9) 



where we have introduced the shorthand notation 



(3.10) 



The expression for TZ(t) can be simplified by introducing the functions 

h(u,t) = C dTe~ iuJT g(r) 
Jo 

k(u,t) = I dr e- iu)T g(r) (3.11) 
Jo 

= iuoh{uo,t) + e- iut g{t). (3.12) 
In terms of these functions, lZ(t) can be written as 

= IhE— (l + 2N(u k ))(\k(u; k ,t)\ 2 + Q 2 \h(u k ,t)\ 2 ). (3.13) 

k 

The expectation value of the number operator ( |3.9| ) in the non-equilibrium density matrix 
has two contributions: one that is completely determined by the initial state of the system 
(proportional to pi ; qi ; a) and the other, determined by the bath and given by 1Z{t). 
Detailed understanding of the particle number relaxation requires the knowledge of the 
dynamical function g(t) which will be studied in the following section. 



B. Calculating g(t) 

Before specifying a choice of the spectral density of the bath J(uj) we can obtain more 
insight by analyzing the real time behavior of g(t) and consequently of < n(t) > in general. 
Having determined the general features of the evolution, we will then specify a particular 
choice of J{w) and provide a detailed numerical study comparing with different approxima- 
tions in a later section. In general the spectral density fullfils 
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1 otherwise 

where ui t h ] oj c are threshold and cutoff frequencies respectively. 

The real time evolution of g(t) can be obtained by taking the Laplace transform of eq. 
(3/7). Solving for the Laplace transform of g(t), namely g(s), one can show that 

s 2 + tuft + S(s) 

with the Laplace transform of the retarded self energy given by 

where we have taken the limit of a continuum distribution of bath oscillators as given by 
The function g(t) is then given by the inverse Laplace transform 

g(t) — ( e st g(s)ds, (3.17) 



27ci Jr 

where T refers to the Bromwich contour running along the imaginary axis to the right of 
all the singularities of g(s) in the complex s-plane. Therefore we need to understand the 
analytic structure of g(s) to obtain the real time dynamics of the particle occupation number. 

From the expression (|3.16|) for the Laplace transform of the retarded self-energy, we find 
that Ss(s) has cuts along the imaginary s-axis for s = iuj ;u>th < M < oj c as can be seen 
from 

£s(s = iuj ± + ) = ± i S/(cj) 

with 

Z R (u)=*V [du>^4 2 (3.18) 

7T J UJ 2 — Uj' 

S/(o;) = 2 sign(cij) J du>'J(u>') u>' 5{oj' 2 — oj 2 ) 

= sign(^) J{\u\) . (3.19) 

It is convenient to introduce a renormalized frequency by performing a subtraction of the 
self-energy. Clearly the subtraction point is arbitrary, and we choose to subtract at s — 0. 
We thus introduce the renormalized frequency as 

2 2 , W„ n\ , .2 ^ f _ J(a>) 



w£ = + £(s = 0) = / rfcj , (3.20) 

71 J OJ 

and the once subtracted self energy is given by 

2 f» , J(oj) s 2 



2 /"OO 

E s (s) = E(s) - S(s = 0) = - / dcj 

7T JO 



OJ S 2 + OJ 2 



10 



Isolated poles of g(s) are at the values s p which satisfy 

4 + uj2 R + ^s(s p ) = 0. 

These are purely imaginary when they are below the threshold frequency of the bath (o;^)(see 
eqn. |3.14|) , corresponding to renormalized exact stable states of the particle-bath interacting 
system. 

If the imaginary part of the pole (in the s- variable) u p is above threshold (uj p > u t h, then 
the pole is in the second (unphysical) Riemann sheet and for weak couplings the spectral 
density S(u), defined below, will feature a Breit-Wigner resonance shape where the width 
of the resonance is related to the imaginary part of the kernel S5 and the peak of the 
resonance is at uj p . The position of these complex poles can be parametrized in terms of real 
and imaginary parts as 

s p = iujp — r. 

These correspond to decaying states and are not eigenstates of the interacting Hamiltonian. 
If the width T << uj p these long-lived resonances are almost energy eigenstates and will be 
identified with the quasiparticles of the interacting system in the next section. 

Depending on the strength of the coupling with the environment, J(uj), and the value of 
ur, the imaginary part of the pole, lj p , can be above or below the threshold, uj th . 
I) Consider first the case in which the pole is above threshold, i.e. u p > uj^. Since there 
are no isolated singularities below threshold, only the cut will contribute to the integral 
( p. 171 ). The Bromwich contour V in the complex s-plane is chosen as the one shown in fig. 
|l](b) where all the singularities of g(s) are to the left of the contour. Evaluating the integral 
along this contour, we obtain 



g{t) = - f U ° dujS{uj) sin(cut), (3.21) 
where the spectral density S(u) is given by 



S(u) = Ej(oj) \g(s = zoo + e)| 2 = — 2 y^L , rv 7 ™ ■ (3-22) 

From the initial condition g(0) = 1 we find the sum rule 

2 r Uc 



dLuS(uj) = l. (3.23) 

For weak coupling, the spectral density can be approximated by a Breit-Wigner resonance 
and asymptotically g(t) is approximately given by pi] 



g[t) ~ Z cos(u p t + a) e ; T 



2u p 



duo 2 



(3.24) 



with a a constant phase-shift [p4]| . We identify this behavior with a typical quasiparticle 
which acquires a width through medium effects and whose residue at the quasiparticle pole, 
i.e. the wave function renormalization constant, is smaller than one as a consequence of 
the overlap between the initial bare particle state and the continuum of the bath. This 
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interpretation will be further clarified when we study the exact normal modes in the next 
section. 

II) Consider next the case in which there is only a single isolated pole below the cut. In 
this case, there are two contributions to the integral (|3.17|) ; the pole contribution and the 
cut contribution. In this case we find 

2 r u <= 

g(t) = Z cos(u p t) H — / cLuuj S(uj) cos(ut) , (3.25) 
where we define the wave function renormalization Z as in fl3T24l) above, 

" ■ (3 ' 26) 

Asymptotically at long time, the cut contribution vanishes with a power law determined by 
the behavior of S(uj) near threshold @], and g(t) oscillates with the pole frequency uj p . Just 
as in the previous case, the bare particle has been dressed by the bath, and to distinguish 
from the bare or quasiparticle we call this state the dressed particle. The position of the 
dressed particle pole has been shifted and its residue is smaller than one as a result of the 
overlap with the continuum of states of the bath. 

From the initial condition g(0) = 1, we derive the important sum rule 

Z + - r duuS(u) = 1. (3.27) 

7T Jui th 

Both cases of the sum rule (|3.23|) and (|3.27|) are a consequence of the canonical commutation 



relations Since the spectral density S{uj) is positive semi-definite, the above sum rule 
determines that Z < 1. 

The expression Q3.25| ) allows us to explore the concept of the dressing time of the particle. 
At long times the contribution to g(t) from the continuum vanishes typically as a power law 
determined by the behavior of the spectral density near threshold [34]] and the contribution 
from the pole dominates the dynamics. This contribution results in a asymptotic oscillatory 
behavior of g(t) with an amplitude determined by the residue Z at the particle pole. The 
formation time can be defined to be the time it takes for the amplitude of g(t) to reach 
its asymptotic value Z (initially g(0) — 1). In the case in which the pole is embedded 
in the continuum (unphysical Riemann sheet) and we deal with quasiparticles, a similar 
concept can be introduced, now being the formation time of the quasiparticle. There are 
now two competing time scales: the formation time scale during which the quasiparticle 
pole dominates the dynamics and the contribution of the continuum becomes subleading, 
and the relaxation time scale which is determined by the imaginary part of the self energy 
at the quasiparticle pole, i.e. the width of the resonance. The time scale of formation of 
the quasiparticle can be defined to be the time it takes until the exponential fall-off of the 
correlation function ensues. 

In this case, the two different time scales can only be resolved if they are widely separated 
which requires that the resonance be very narrow and the exponential relaxation associated 
with the decay of the quasiparticle allows many oscillations to occur. This condition can be 
quantified as V /u p « 1 which requires a weak coupling to the bath. We will explore these 
situations numerically in a later section where a particular density of states of the bath will 
be proposed. 
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C. Fluctuation-Dissipation 



The main advantage of studying the time evolution of the reduced density matrix is 
that it allows to establish a direct relationship between the relaxation of the occupation 
number of the "system" and the fluctuation dissipation theorem. The connection between 
the fluctuation-dissipation and the Boltzmann equation has been investigated recently in 



the semi-classical regime [13 



in terms 



This relationship is established by re- writing the path integral (P.I.) in eq. 
of the Wigner coordinates which can be cast in the following probabilistic form p7| , |28| , |30 

P.I. = f VxVrVtPl^e^??^^ 

S eff (x, r, = J* dt'r(t') j- (x(t') + uj 2 x(t')) + J dt"Z(t' - t")x{t") + £(*') 

P\i] = exp {-i JdtJ dt>at)K-\t - *')£(*')} • (3-28) 

The path integral over the relative coordinate, r(t), leads to a non-Markovian Langevin 
equation for the center of mass coordinate, x(t), in the presence of a stochastic Gaussian 
(but colored) noise term £(£) [E7,28,30|. The noise correlation function is determined by 



K{t — t') given by eq. (|3.6|) . 

The fluctuation-dissipation relation is established in the following manner pi |. In the 
limit of a continuum distribution of the bath oscillators we find the time Fourier transform 
of the retarded self-energy £(t), eq. fl3.8|), to be given by the analytic continuation of the 
Laplace transform (|3.16Q s — > uj — ie, i.e. 



Then we find (for u > 0) 



ie 



-- I ,U 
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J J{J) 



\{uj') 2 -(lu- ie) 2 ] 



Im 



J(u) 



and the Fourier transform in time of the kernel K(t) (|3.6|) is given by 



K(u) = — Im 



coth 



0u 



(3.29) 



(3.30) 



(3.31) 



This is the usual fluctuation-dissipation relation [^TJ. Finally we obtain the bath contribution 
to the non-equilibrium occupation number eqn.( [3.9|) , which is determined by TZ(t) given by 
eqn.( ^.13] ), in a form that displays clearly its relationship to the fluctuation dissipation 
relation 



K{t) = ^ / dwK{u) ( \k{u,t)\ 2 + Q 2 \h{u,t)\ 2 



(3.32) 



where K(u) is the power spectrum of the bath. This expression makes explicit the stochastic 
nature of thermalization and establishes a direct relationship with the fluctuation dissipation 
theorem. 
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IV. THE HEISENBERG OPERATORS 



The above results can be understood in an alternative manner by obtaining the real 
time evolution of the Heisenberg picture operators, from which the expectation value of the 
number operator can be obtained by providing an initial density matrix. 

The equation of motion of the Heisenberg operator q(t) is given by 



q(t) + J*q{t) - £ 3 f dt > sinMt-t')] = - J2 C kQ k 0) (t), C 4 - 1 ) 



where Q k (i) satisfies the homogeneous equation 



(o) 



Qr it) + uiQt\t) 







(4.2) 



Eq. (|4.1| ) can be solved using Laplace transform, and the operator solution with the initial 
condition q(t = 0) = q(0) ; q(t = 0) = p(0) is found to be given by 



q(t) = p(0)g(t) + q(0)g(t) - E C k f dr (t - r) g(n 

, Jo 



(4.3) 



where g(t) is the same function which was defined in the previous section. 

Since the initial density matrix describes a thermal distribution for the quanta of a 
harmonic oscillator of reference frequency Q, it is convenient to write the initial position 
and momentum operators in terms of the creation and annihilation operators of a quanta of 
frequency Q as 



q(0) 



1 



p(0) 



^2 



Also, it is convenient to write Q^\t) in terms of the creation and annihilation operators of 
a quanta of frequency ojk as 



Gathering all terms, q(t) and p(t) become 
1 



a k e 



t „iui k t 



al e 



pit) 



2tt 
1 



Ck 



b (g(t) - iSlg(t)) + 6 f (g(t) + i{lg(t))\ - E , 
b (g(t) - Xlg(t)) + 6 f (g(t) + itlg{t))] - E ( ' 



i{ e iuJkt h(u k ,t) + h.c. 



ale iuJkt k(to k ,t) + h.c. , (4.5) 



(4.4) 



where h(ui k , t) ; k(cv k , t) are defined in eq. ( |3.12| ). The expression ( |4.5[ ) reveals that the particle 
operators create states with overlap with bath continuum. 

The expectation value of the occupation number operator n(t) in ( |2.2p can be evaluated 
using an initial density matrix which is diagonal in the basis of the number operators for 
system and bath. Assuming a continuum spectrum of the bath oscillators, using ( p.l|) and 
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considering for simplicity the case of vanishing expectation values of q(0) ; p(0) in the initial 
density matrix, we find 



<«(*)> = - ^ + {m + (t) + ny (*)) 

+ ^JduJ(u)(l+2N(oj))(\k(uj,t)\ 2 + Q 2 \h(oj,t)\ 2 ). (4.6) 

Setting qi = ; pi = in the result (|3.9|) we find that ( [4.6| ) reduces to the expression 
obtained by the time evolution of the density matrix ( |3.9|) and the last term is identified 
with K(t). 

The operator method allows to compute any correlation function of operators in the 
initial density matrix at arbitrary times, whereas the time evolution of the density matrix 
would require the introduction of external sources and taking functional derivatives with 
respect to those to obtain unequal time correlation functions. 



V. ASYMPTOTIC BEHAVIOR OF THE OCCUPATION NUMBER 

The asymptotic behavior of (n(t)) is completely determined by the long time dynamics 
of g(t). We have shown that g(t) vanishes asymptotically for poles in the continuum while 
the contribution from the isolated pole dominates for the case in which the pole is below 
threshold. We will consider each individual case in detail. 



A. Poles in the continuum (lo p > uith) 

In this case the function g(t) vanishes exponentially at asymptotically long times (ft.24|) 
and the asymptotic behavior of the particle occupation number is given by 

(n(oo)) = -i + fc(oo) (5.1) 

with 

H{oo) = —— r° du \l + 2N(u)} S(u) (Q 2 + uj 2 ) , (5.2) 



where we used eqs.( |3.19| ) and ( |3.22j ), recognized the Laplace transform of g(t) in the long 
time limit for eq.( p,13| ) (using the vanishing of g(t) at long times) 

\h(u, oo)| 2 = \g(s = iuo + e)| 2 
\k(u, oo)| 2 = u 2 \h(u, oo)| 2 . 

It is clear that the asymptotic value of (n(oo)) is different from the equilibrium occupation 
number of the bath N(u p ). 

Suppose that the spectral density S(lj) can be approximated by a narrow Breit-Wigner 
resonance with 
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S ( U ) = 7TT T~. \o i -n? ^ 7T— ( 5 - 3 ) 



7rZ 



2u; p (a; — u p ) 2 + T 2 2lu p 
where 

ZS/(u p ) 



(5.4) 



as would be the case for weak coupling. Then the asymptotic occupation number becomes 

1 



(n(oo)) = Z 



+ - 



2 ' (5 ' 5) 



which is different from the equilibrium value of the bath. 

We now see that choosing the reference frequency fl = u p , and introducing a quasiparticle 
number operator as 



we see that 



n q uadt) = ^[p 2 (t)+uj 2 p q 2 (t) -— (5.6) 



1 

2Z 



N{u p ) + \{1-^) (5.7) 

thus the quasiparticle reaches an asymptotic distribution function of almost thermal form 
with corrections arising from the wave function renormalization. The factor Z in the defi- 
nition of the quasiparticle number operator reflects the fact that the quasiparticle pole has 
strength Z rather than one. In the following section it will become clear that the result 
( |5.7| ) is a consequence of the complete thermalization of the quasiparticle with the bath and 
that the occupation number of quasiparticles becomes the one predicted by the canonical 
ensemble as follows from the discussion leading to the eqs.( |6.9| , |6.10| ) below. This expression, 



thus reveals the importance of counting the quasiparticles instead of the bare particles. Even 
in the weak coupling limit the distribution of bare particles is not thermal whereas the true 
quasiparticle distribution departs perturbatively from a Bose Einstein distribution at the 
temperature of the bath. 

The asymptotic value of the distribution is approached exponentially. The thermalization 
time scale is given by = 1/21" since it is determined by g 2 (t) which is the dependence of 
the occupation number on the function that determines the real time evolution either of the 
density matrix or of the Heisenberg operators. 

Even when the occupation number is defined in terms of the true 'in medium' pole, there 
will be departures from the Bose-Einstein distribution for non-negligible width T and when 
the strength of the pole Z is substantially smaller than one. These corrections will arise 
in the case of broad resonances and may lead to large departures from the Bose-Einstein 
distribution. This situation will be explored numerically later. 

In the case of a wide resonance, the product N(cu) S(cu) is sensitive to the width of the 
resonance. For bath temperature T « u) t h the Bose-Einstein distribution will only probe 
the tail of the broad spectral density closer to threshold and the product is only sensitive to 
the threshold behavior of S(u) 
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In particular if near threshold S(u) ~ (u — uJth) a then for temperatures T « ujth the 
temperature dependence of the equilibrium abundance of unstable particles in the bath is 
approximately given by 



n(T; t = oo) — n(0; t = oo) « e r 



-1Q + 1 



which reveals threshold corrections to the Boltzmann exponential suppression. This result 
has been anticipated in within a different context. 



In the opposite limit when T » uj p the product is sensitive to the width of the resonance 
and the details of the spectral density. Thus in the case of a broad resonance the departures 
from a Bose-Einstein distribution function for the quasiparticles will be important. Clearly 
this is the regime in which a Boltzmann approximation could be unreliable. 

These corrections originate in off-shell effects that will depend on the particular spectral 
density of the bath and the coupling between the particle of the bath. We will quantify the 
corrections for a particular choice of J(u) in a following section. 

Moreover, the asymptotic value of the particle number does not depend on the initial con- 
dition of the particle; e.g. initial expectation values of position and momentum, temperature 
or occupation number. 



B. Isolated poles (lo p < uj t h) 

In this case the asymptotic time dependence of the function g(t) is completely determined 
by the isolated pole below the bath continuum and the function "rings" with this frequency 
and with asymptotic amplitude determined by the wave function renormalization Z given 
by (|3.26|) . The asymptotic behavior of the particle occupation number defined at a reference 
frequency Q is now given by 

r72 "~ 2 ' 'p^ + ^ + ^a + ^qf 

PiQiZ 2 ( Vt uj„\ , , Z 2 cos 2 (uj ri t) ( pf _ _ 9 \ . 
+ -^{- r -t) Shl(2 ^» + 2^ (t + 2!V + n "'j • (5 ' 8) 

where 72. (oo) is the limit value of TZ(t). For Q = u p , i.e. the position of the dressed particle 
pole, the asymptotic value of the occupation number obtains the simple form 

(n(oo)) = -- + K(oo) + Z 2 

The last term can be identified as the contribution from the expectation values of 
p(0) ; q(Q) (pi ; qi respectively) in the initial density matrix. 

Unlike the case in which the pole is in the continuum, the asymptotic value of the 
particle occupation does depend on how the particle was prepared initially since expression 
Q5.9D depends on p^, q^ and n(0). 

In this case, 72. (oo) has contributions from both the continuum cut and the isolated pole 
below the continuum. 

In order to compare the results to those obtained from an approximate quantum kinetic 
equation obtained via a perturbative expansion in the next section, it is useful to obtain an 



M°o)> = -- + 72(oo) + \ 2 pJ 



n(0) + 



+ 



2fi 



(5.9) 
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expression for TZ(oo) up to first order in J{uS). The expression for TZ(oo) ( |5.2| ) is proportional 
to the spectral density 5* (a;) given by ( p.22| ) . When the pole is below the continuum, the con- 
tribution from the cut is proportional to J{uj) and perturbatively small when J(uj) is small. 
Furthermore the continuum contribution dephases rapidly at long times, and asymptotically 
the relevant contribution to g(t) arises from the isolated pole. After some straightforward 
algebra we find for Q = u p that at long times 

i,/ si9 ^i,/ m9 ^9 f 1 ~ cos(u;+t) 1 — cos(u^t) 1 /_ „„\ 

k(u,t)\ 2 + n\h(u,t)\ 2 = Z 2 { 9^-^ + t — -t 5 - 10 

[ v+ u- J 

with 



u± = Q ± uj 

and to lowest order in J(u), the asymptotic contribution TZ(oo) is given by 



(5.11) 



TZ(oo) 



2ntt 



du J(uo) [1 + 2N(u)] 



1 



1 



o(j 2 ) 



which for easier comparison with the results from kinetics, can be written in the following 
form 



1 + N(u) N(u) 



G(J 2 ) 



where the term (1 — Z 2 ) m 2(1 — Z) and we have used the sum rule ( |3.27| ) to lowest order. 
Setting pi = qi = in eq.( |5.9|) , the asymptotic occupation number becomes 



(n(oo)) 



n(0) + 4^ / ^^(^) 

7TSi j 



1 + A^(cj) iV(a;) 



^4 



UJZ 



+ J 



(5.12) 



We have purposedly kept Z in the above expression to compare it to the results from the 
quantum kinetics approximation to be obtained later. 

Clearly this result depends on the initial distribution of the particle and the details of 
the spectral density of the bath, leading to the conclusion that in the case in which the 
particle pole is real (below threshold), the particle does not thermalize with the bath. 



VI. COLLECTIVE NORMAL MODES AND QUASIPARTICLES 

In a many body problem, the poles of the exact two particle Green's functions are 
identified with the collective modes. In general the poles are complex resulting in the 
damping of the collective excitations. We can make contact with this many body concept 
by studying the normal modes of the total Hamiltonian for the particle-bath system under 
consideration. 

Since the Hamiltonian is quadratic, it can be diagonalized by a canonical transformation 
in terms of the normal modes. In order to establish a correspondence with the continuum 
distribution of bath oscillators it is convenient to write the Hamiltonian in the continuum 
form 
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H 



£(p 2 + " V) 



duj 



P 2 {uj) + uj 2 Q 2 {uj) +q duo C{uj) Q(uj 
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2uj 



The Hamiltonian of this rather simple model can be diagonalized by finding the normal 
modes. Let us write the linear change coordinates and momenta (canonical transformation) 
to the normal modes as 



q = S\ a(X)Q(X) ; p = S x a(\)T(\) (6.1) 
Q(u) = S x P(u,\)Q(\) ; P(u)=S x P(u,\)P(\), (6.2) 

where the symbol S x stands for the sum over the discrete and integral over the continuum 
normal mode eigenvalues A that render the Hamiltonian in diagonal form 

1 



H 



-S x 



V 2 (X) + X 2 Q 2 (X) 



The vectors V(X) = (a(A), (3(ui, A)) obey the normal mode eigenvalue equation which in 
components reads 



uj 2 a{X) + / duj C{uj) (3(uj, A) = A 2 a(A) 



C(cu) a(A) + uj 2 P(cu, A) = A 2 0(u, A) 



(6.3) 
(6.4) 



and the X's are the exact eigenenergies of the Hamiltonian. 

Solving for (3{u, A) in terms of a(A) in eq. ( |6.4j) and inserting the solution back into ( |6.3|) 
we find the solution for the coefficients and the secular equation for the eigenvalues to be 
given by 

C(u) a{\) 



(A 



ie) 2 — uJ 2 



+ B 5(X-lu) 



A - uj Q 

7T 



duj— 



UJ J(lj) 



a{X) = B C(A) 



where we used 'retarded' boundary conditions (with the ie prescription) to establish contact 
with the previous results, and B is determined by normalizing the eigenstates. 

There are two distinct possibilities: I) an isolated pole below the continuum threshold 
of the bath corresponding to a dressed stable particle, II) a continuum of states and a 
quasiparticle pole in the unphysical Riemann sheet (resonance). 

I) Isolated poles: The condition for isolated poles below the bath continuum requires 
setting B = since the spectrum of the bath has no support below threshold. The position 
of the pole is found from the secular equation 



2 2 



7T Jw t 



du) 



u J{uj) 



. 



~-Jt 



UJ^ 



This expression is identified as the condition for isolated poles in the Laplace transform g(s) 
(see eq. (|3.15p ) for s 
find 

a(uj p ) = \fz 



iujp. The value of a(oj p ) is determined from normalization and we 
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with Z the wave function renormalization given by eqs. (|3.18|) and (|3.26|) . Normalization of 



the vectors is equivalent to the sum rule ( 3.27 ). 
II) Continuum states: 

For the continuum states we take B = 1 so that Q(A) — > Q(u) when C(u) — > and we 
find the coefficients 

a(A) = rr-r- (6.5) 

^A) = 5(X-u) + { C^l 2 . (6.6) 

In this case the normalization results in the sum rule given by eq.( [3.23| ). 

Because of our choice of boundary conditions, the coefficients are complex and the re- 
sulting new coordinates are not Hermitian. This can be remedied by absorbing the phases 
by a trivial canonical transformation and defining the coefficients a (A) ; /3(A) in terms of 
their absolute values and the Q(A) ; V(X) to be real. This phase carries the information 
of the boundary conditions (the ie prescription) and since it is removed by a canonical 
transformation the results are independent of these. 

Let us consider the case of an isolated pole below the threshold of the bath continuum 
at A = uj p . This state is the one that evolves from the bare particle degree of freedom upon 
adiabatically switching-on the system-bath couplings and is identified with the position 
of the isolated pole in the Laplace transform of the function g(t) given by (|3.15| ). 

Separating the contribution from the isolated pole we write 

q(t) = VZQo(t) + Q con t(*) 
p(t) = VzV (t)+V cont (t) , 

where the operators <2 con t ; ^cont crea te excitations in the continuum of the bath out of 
the exact ground state. Writing Q Q ; Vq in terms of creation and annihilation operators of 
the exact eigenstates, we see that asymptotically long times the operators q(t) ; p(t) create 
an exact one dressed particle state out of the exact vacuum. In the limit of asymptotically 
long times and invoking the Riemann-Lebesgue lemma 

ff(*)|0> - ^I=e^|lp>, 
V AuJ p 

where |0 > is the exact ground state and the contribution from the continuum states averages 
to zero at long times by the dephasing between modes. 
The operator 

A\{t) = ^L=[u p q{t)+ip{t)} (6.7) 
2ujpZ 

asymptotically at long times creates a dressed particle state with unit residue out of the exact 
vacuum. At any finite time the state created by this operator is not an eigenstate of the full 



Hamiltonian but has overlap with states in the continuum [35]. We associate the operator 
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( p.7| ) with dressed particles in the case of isolated poles or quasiparticles for resonances, in 
contrast to the normal (collective) modes of the system that are exact eigenstates. 

Although a priori one would be tempted to define the dressed particle as the normal 
mode of frequency uo p associated with the creation and annihilation operators obtained from 
the normal mode described by Qq ; Vq, these are of little use: these operators represent 
linear combinations of the particle and the degrees of freedom of the bath. Obviously the 
number operator associated with this normal mode is constant in time. The interpolating 
operator ( |6.7| ) is the natural candidate for counting quasiparticles . 



In an experimental situation such as for example an electron in a metal, one would 
like to write down an evolution equation for the distribution function that describes the 
particle dressed by the medium. The interpretation of the quasiparticle creation operator is 
consistent with this physical situation since the added particle will move in the bath being 
dressed by the interaction with the medium, the resulting quasiparticle will have a new 
dispersion relation (given here by u p ) and in general a width, and the probability associated 
with this quasiparticle pole will be reduced by the overlap with the states of the bath. This 
quasiparticle is not a stationary state because it overlaps with the collective modes and its 
time evolution involves dephasing. 

In the case in which the pole at u p has a value larger than the threshold for the bath 
oscillators, it has moved into the second (unphysical) Riemann sheet upon adiabatically 
switching-on the interaction and is no longer part of the eigenspectrum of the Hamiltonian. 
In this case it has become an unstable state and uj p will have a small negative imaginary 
part given by V [see eq. (|3.24j )] . In this case the overlap with the continuum results in an 



almost exponential decay of the quasiparticle distribution after the short formation time of 
the quasiparticle. 

We then see that the interpolating number operator 



2u p Z 



h q uasi{t) = tt^ P\t) + u 2 p q 2 (t) - (6.8) 



2Z 



can be interpreted as either the dressed particle distribution function in the case of an isolated 
pole below the continuum of the bath or the quasiparticle distribution function in the case 
of a resonance. Besides setting the reference frequency A = io p in eq.( |2.2| ) the wavefunction 
renormalization factor Z accounts for the strength of the particle or quasiparticle pole. The 
importance of wave function renormalization has been highlighted within the context of high 



field transport in semiconductors |36] 
Interpretation of results: 

This analysis in terms of normal modes reveals several features of the exact solutions 
obtained in the previous sections. 

• Thermalization of quasiparticles (resonances): in the case in which the quasi- 
particle pole is above threshold, the asymptotic value of the quasiparticle distribution 
given by eqs.( |5.1D -( |5.2|) is a consequence of thermalization. Indeed by using the ex- 
pansion of q ; p in terms of the normal mode coordinates and momenta given by 
eqs.(|6~T| ) -(|6.2| ) with the coefficients |a(A)| and real V(\); ; Q(A) it is straightforward 



to prove that 

1 

Z 



< n quasi (oo) >= Tr 



f^quasi (0) 



-I3H 



K(oo) - \ 



(6.9) 
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with n quas i(0) the quasiparticle number operator Q6.8Q at the initial time t = 0. This 
is a remarkable result: the density matrix, which initially was of a factorized form 
for particle and bath at different temperatures has evolved in time to the equilibrium 
density matrix for the total system at the temperature of the bath. However the 
distribution of quasiparticles is not given by the Bose-Einstein form. Furthermore, 
the contribution to n quas i(oo) that does not vanish as T —>■ can be interpreted as 
a zero point contribution from the resonance. In the case in which the quasiparticle 
becomes a narrow resonance we see from eqs.( |5.5|) and ( |6.8| ) that 

id-I 

2 V Z 



< n quasi (oo) >= N(u p ) + -(1 - -) (6.10) 



and the number of quasiparticles departs from a Bose-Einstein distribution at the 
temperature of the bath with the departure determined by the off-shell effects that 
result in Z ^ 1 through the sum rules. The last term, identified above with the zero 
point contribution is interpreted as the normalization borrowed from the continuum 
by the quasiparticle. Although in this simple case Z does not depend on temperature 
and the last term in (|6.10|) can be subtracted out as a redefinition of the quasiparticle 
vacuum, in a general quantum many body theory, the wave function renormalization 
will be medium dependent and such subtraction would be unjustified. 

Non-Thermalization of stable particles: In the case in which the particle pole 
is below threshold the asymptotic oscillations in the expression ( |5.8|) for Q ^ uo p are 
a consequence of the interference between the state of arbitrary frequency Q and the 
normal mode with frequency u p . These oscillations disappear when the reference 
frequency (Q) is chosen to be the normal mode pole frequency (uj p ) which is also 
the particle frequency, this fact has already been noticed within a different context 
[Q. The factor Z 2 in eq.( |5.9|) has the following origin: asymptotically at long times 
q(t) -> v^Qo ; pit) -> y/ZV in the sense of matrix elements. But the Qo ; Vq create 
particle states out of bare states with amplitude \f~Z, therefore one of the factors Z in 
eq.( |5.9| ) arises from the asymptotic (weak) limit on the operators, and another factor 
Z arises because the calculation of eq. (|5.9|) was performed in terms of the bare states 
overlap with the particle states given by the wave function renormalization. Using the 



expansion in terms of normal modes we find that < n quasi (oo) > given by eq.( |6.8|) does 
not coincide with 



Tr 



unlike the previous case of a resonance. 



VII. KINETICS 

Having provided an analysis of the exact evolution of the distribution function and dis- 
tinguished between that of renormalized, stable particles and quasiparticles (resonances), 
we now proceed to obtain kinetic equations in several stages of approximation to compare 
with the exact results. Kinetic equations are obtained by truncating the hierarchy of equa- 
tions of motion for the higher order correlation functions under certain assumptions. The 
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typical assumptions are those of slow relaxation as compared to the microscopic time and 
length scales and rely on a separation of scales. To warrant this separation between scales 
clearly a perturbative parameter must be invoked and the resulting kinetic equations pro- 
vide a resummation of the perturbative expansion. Different type of approximations result 
in different resummation schemes. 



A. Boltzmann Equation 

The simplest kinetic equation to describe the approach to equilibrium is the Boltzmann 
equation, which is obtained by writing a gain minus loss balance equation in which energy 
is conserved and using Fermi's Golden Rule. Writing q(t) and Qkif) m terms of creation 
and annihilation operators, we identify the only terms that can contribute by energy con- 
servation: the 'gain' term in which one quanta of the system's oscillator is created and one 
quanta of the oscillator with label k is annihilated, minus the 'loss' term in which a quanta 
of the system's oscillator is annihilated and a quanta of the oscillator of label k in the bath 
is created. The first term has probability given by 

C 2 

gain = — -*— (l+ n )JV fc . 
4cu uj k 

The second term has a probability 

loss= - k (1 + N k ) n. 
Thus using Fermi's Golden Rule we find the usual Boltzmann equation 

(h B (t)) = (2tt) ]T [(1 + n(t))N k - (1 + N k )n(t)} 5{u h - u ) -> 



duo [(1 + n(t))N(u) - (1 + N(uj))n(t)} 5(uj - w ), (7.1) 

where we have assumed that the bath remains in thermal equilibrium with constant distri- 
bution functions. The solution is clearly 

<n fl (f)> = NM + [n(0) - N(u )} e"*™ ;T B = ^ = (7.2) 

where we recognize the lowest order (Born approximation) to the decay rate which is given 
by eqs.( |3.19| ) and (|3.24j) . Obviously the Boltzmann equation predicts no relaxation in the 



case in which the pole is below the continuum threshold since in this case J(ujq) = 0. 

Even when the bare frequency is in the continuum of the bath, the Boltzmann approxi- 
mation predicts no relaxation if n(0) = N(uj ) as the gain and loss processes exactly balance. 

As we will see explicitly numerically below the exact solution shows a non-trivial time 
dependence in this case because the bare particle is dressed by the medium and the asymp- 
totic equilibrium distribution function reveals off-shell effects as discussed in the previous 
section. 
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B. Quantum Kinetic Equation: 



The quantum kinetic equation is obtained by taking the expectation value of the number 
operator using the Heisenberg equations of motion and truncating the exact equations of 
motion within a particular approximation. 

Since we want to obtain the kinetic equation for the relaxation of the distribution function 
of particles with frequency Q (for quasiparticles this is the pole frequency of the propagator, 
for bare particles it is simply u ) it is convenient to write the total Hamiltonian in terms of 
this frequency adding a counterterm of the form 

x, ,2 

H ct = ^q\t) ; 5u 2 = u 2 -Q 2 . 

As usual the counterterm is chosen appropriately in perturbation theory to cancel the con- 
tributions recognized as those arising from a shift in the frequency. 

Taking the derivative of eq.( [2.2|) and using the equations of motion we obtain 



*(*) = -ji {ec* Qmm + ^ im m+m m}] ■ (7.3) 

The expectation value of the time derivative of the occupation number is calculated by 
multiplying eq. (|7.3| ) by p(0) and taking the trace 

..,2 1 

(7-4) 

t'=t 



m = -~ fee* (q + (t') Q k it)) + 5 4- mm) + mm) 



q + (t')Qk(t)) =Tr [q(t')p(0)Q k (t)] . 



Qdt 

where 



We need to evaluate the non-equilibrium matrix element (^q + (£') (t)) . This can be 
achieved by treating the interaction term in perturbation theory. The zeroth order term 
in the perturbative series does not contribute because the initial density matrix commutes 
with the number operator at the initial time. 

A simple diagrammatic analysis of the perturbative series reveals that the kinetic equa- 
tion can be written exactly as 



n(f)> = -^E^|; { \fdt"(E<(t,t")g>(t",t') - £>(M")£<(£V)) 

h 



+ 



X, ,2 



where Q <,y are the exact Green's functions for the system, defined similarly to those of the 
bath eq.( |3.1|) and S fe ' are the irreducible self-energy components, again defined similarly 
to eq.fni). 

To first order in the interaction we use the free-field propagators and the lowest order con- 
tribution to the self-energy. It is straightforward to show that the counterterm contribution 
vanishes to this order and eq.( [7.4[ ) becomes 
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(«(*)> =h^° l M [J™ dt" ((q + (t') q +(f))(Q^t)QUO] 

- (q + (t')q-(t"))(Q k (t)Q k (t")))] t , =t ■ (7.5) 
Substituting the non-equilibrium Green's functions from eq.( |3.1| ) in the right hand side of 



eq. (17.50 ) taking the derivative with respect to t' and arranging terms, we obtain 

(n(t)) = — f dt' I duoJ(uo) {[1 + ti(0) + N(u)\ cos \(Q +u)(t- t')] 
nil Jo J 

+ [N(u) - n(0)] cos [(Q -u)(t- t')} } , (7.6) 

where n(0) is the distribution of quanta for the particle at the initial time and N(u) are 
the Bose-Einstein distributions of the bath which will be taken to be constant and given by 
eq.fl3j). 

We now propose a scheme that provides a resummation of the perturbative series by 
replacing the initial distribution n(0) by self-consistently updating the distribution inside 
the integral in eq.( |7.6| ) by replacing n(0) — * n(t'). It will be shown explicitly below that 
this prescription leads to a Dyson summation of particular Feynman diagrams and the case 
where n is constant is understood as the lowest order term in this expansion. Within non- 
relativistic many-body quantum kinetics, this approximation is known as the generalized 
Kadanoff-Baym ansatz ||17|| - ||22|| . The validity of this approximation in the weak coupling 



limit is confirmed by comparing the resulting evolution of the distribution function to the 
exact result obtained in the previous sections as it will be seen below in detail. 
The quantum kinetic equation is then given by 

(n qk (t)) = — J dt' J duo J{uj) {[1 + n(t') + N(u))\ cos [(CI + u)(t- t')} 

+ [N(u) - n(t')\ cos [(tt -u)(t- t')} } . (7.7) 

The resulting linear kinetic equation, eq. (|7.7| ), can now be solved via Laplace transforms. 
The Laplace transform of < n q k(t) > is given by 



s , N(lj) 

.■2 4-/. i 2 X 



n q k(s) = f ^ -, (7.8) 



s 

T 



where n(0) is the initial occupation number of the particle and u± are given by eq.( ^.ll| ). 
The dynamics of the occupation number of the particle is obtained by taking the inverse 
Laplace transform along the Bromwich contour. The analytic structure of n q k(s) consists 
of cuts along the imaginary axis in the s-plane and poles. For the pole contributions, we 
distinguish two cases : 

Case I : Poles in the continuum. In this case there are two poles: 1) a pole where the 
denominator of eq.([7.8|) vanishes, i.e. 



Sp 
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For weak coupling, one can solve for the pole, s p , in perturbation theory and one can show 
that the pole is given by (up to first order in J (a;)) 



s P n 



where we used the identity 



! h S?T^4 =7r ^ ±) ' (7 ' 9) 

and T is given by eq. ( |5.4| ). The contribution from this pole vanishes exponentially for long 
times. 2) There is a second pole at s = and the residue of this pole, using eq. ( [Hp , is 
N(Q). The average occupation number is then given by the contribution of the two poles 
and the cut 

(n(t)) = N(Cl) + (residue at s p ) e~ 2Vt + contribution from the cut . 



Asymptotically the contribution of the cut falls off as a power law Therefore the 

contribution from the last two terms vanish and the particle occupation number approaches 
the equilibrium occupation of the bath with frequency Q. Comparing the above result for 
Q = ojp with the one obtained exactly in the small coupling regime, eq. (|5.5| ), we see that 
they differ by a factor of order J(uj) which can be compensated for by considering higher 
orders in deriving the kinetic equation. Thus we see that for weak coupling (u p « ui ), the 
solution of the quantum kinetic equation approaches that of the Boltzmann approximation 
given by eq.( [7.2|) . 

Case II : Poles below the continuum. Since J{ui p ) vanishes for poles below the continuum, 
there is only one pole at s = 0. The average occupation number is given by the sum of the 
residue of the pole and the cut contribution. At long times the cut contribution vanishes at 
least as a power law [33| and the asymptotic average occupation number is given by 



n(o) + £ r t ; duj jh { 


1+iVM , N(w)\ 






\ 



(n qk (oo)) = 



The denominator of the above equation can be simplified considerably becoming simply Z 2 
to this order. The above equation is now written as 



(n qk (oo)) 



1 p>c ( 1 + N(u) N(lu 



n(0) + duo J(u) < 2^ + 




Comparing the above result with the one obtained exactly in the small coupling regime, 
eq. (|5.12|) , we see that the two results coincide. 

Obviously this quantum kinetic equation includes contributions from intermediate states 
that do not conserve energy and therefore provide off-shell corrections. For the case in 
which the quasiparticle pole is in the continuum, we see that asymptotically at long times 
the distribution becomes similar to that obtained in the Boltzmann approximation with the 
same relaxation rate. However at early times the solution of the quantum kinetic equation 
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differs appreciably from the Boltzmann solution in that the relaxation rate vanishes at the 
initial time, whereas it is a constant for Boltzmann. The vanishing of the relaxation rate 
at the initial time is a consequence of the fact that the initial density matrix is diagonal in 
the number representation, thus whereas the quantum kinetic equation describes correctly 
the initial evolution, the Boltzmann equation has coarse grained over these time scales and 
misses the early time behavior. 

This is important experimentally if the resolution in time of the measurement allows 
to study time scales that reveal features of the initial preparation. Such is the case in 
femtosecond resolved studies of relaxation of hot carriers as described in the introduction. 



C. Markovian approximation: 

If the particle occupation number varies on time scales larger than the memory of the 
kernel in the kinetic equation, a Markovian approximation may be reasonable. In this 
approximation, the particle occupation number n(t') in eq. (|7.7|) is replaced by n(t) and 
taken outside the integral. This approximation would be justified in a weak coupling limit, 
in this case when the spectral density of the bath J(uS) includes a small coupling (as it will be 
specified in the next section) 77. The rational behind this approximation is the realization of 
multi-time scales: a microscopic or short time scale given by t ~ V&'p and another relaxation 
or long time scale t\ ~ rjt. 

Thus in the Markovian approximation, eq.( [7.7[ ) becomes 



1 I A it \ \ (l + ^H)sinKt) 

7T\l J UJi 



+ 



N(u) sin(cu_t) 



+ ^ fduJ{u) | Sin( " +t) - Sin( "- t} I . (7.10) 
nil J [ u + ijj- J 

A computational advantage of this equation is that it provides a local update procedure. 
A connection with the Boltzmann approximation is made with a second stage of approxi- 
mation, known in the Boltzmann literature as the 'completed collision approximation' and 
consists in taking the limit t — > 00 in the arguments of the sine functions in eq.( [7.10| ). This 
approximation enforeces strict energy conservation. This can be understood by using the 
limiting distribution 

mn[w±t] 

hmt^oo = Tid{uj±) 

u± 

which is used in the derivation of Fermi's Golden Rule. Noticing that only u;_ could vanish 
leading to the Boltzmann expression 

(n(f)> = ^ [N(Q) - n(t)] 

which coincides with the Boltzmann equation obtained above |eq. (|7.1|) l within the limit of 
validity of the perturbative expansion, since perturbatively Q ~ c^o if ^ is taken as the 
position of the quasiparticle pole. 
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VIII. NUMERICAL ANALYSIS 



In order to compare the particle number relaxation n(t) between the exact results eg. ( |3.9| ) 
and the various approximations to the kinetic description eq.( |7.7p , Boltzmann and Marko- 
vian, we have solved numerically for a particular choice of the spectral density of the bath. 

We will choose the following model for J(uj) 

J(u) = r] (lj - u t h) d(u - u th ) 9(lu c - u). (8.1) 

This is a generalization of the Ohmic bath in which J{oS) vanishes for frequencies below a 
threshold frequency u t h and above a cutoff frequency uj c , and r\ is a coupling parameter. 
This is the simplest spectral density of the bath that allows us to model important features 
of relevant microscopic models and illuminates the main aspects of relaxational dynamics. 
This form of the spectral density for the bath has been motivated by previous studies of 



decoherence and dissipation in similar model theories ||26|| - [p2]1 . It is the simplest realization 
that allows us to vary parameters and investigate the different regimes for the phenomena 
discussed in the previous section. By varying the coupling r\ and the value of the bare (or 
renormalized) frequency we can test the different scenarios. 

For the case of the quasiparticle pole embedded in the continuum the dimensionless 
parameter that determines the separation of time scales is given for the spectral density 
©by 

r ~ V ,. 



UJth 



When this ratio is << 1 the resonance is rather narrow and there are many oscillations before 
the decay, the time scales are widely separated. In the other limit when this ratio ~ 1 the 
particle is strongly coupled to the bath, resulting in a wide resonance and a potential for 
large off-shell effects including effects related to the proximity of the peak of the resonance 
to the threshold. 

The dynamical function g(t) satisfies eq. ( |3.7| ). In terms of the renormalized frequency 
uir given by eq.( |3.20| ), g(t) can be shown to satisfy the following equation 

g(t)+u 2 R g(t) + - f dt' f dw^-cos[u{t-t')]g{t') = 0; g{Q) = ; g(0) = 1. 
n Jo J to 

We scale our results to an arbitrary unit of frequency and refer all dimensionful quantities 
to this unit since the important physical quantities are dimensionless ratios (such as uj/T 
etc). 

Now we study different scenarios in detail. 

Figure |2] shows the case for which the dressed particle pole is below threshold. In this 
case the Boltzmann equation predicts that no relaxation occurs because the imaginary part 
of the self-energy evaluated on shell (damping rate) vanishes. The exact solution, and 
the quantum kinetic approximation along with the Markovian limit all predict non-trivial 
relaxation, and for this weak coupling case all agree to within few percent. Obviously in 
this case the relaxation is solely due to off-shell effects since the dissipative effects associated 
with processes that conserve energy (on-shell) vanish. The insert of the figure shows the 
dynamics of dressing of the particle and the time scales predicted by the exact result are 
well reproduced by both the quantum kinetic equation and its Markovian approximation. 
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In contrast, fig|| shows the case for which the quasiparticle pole is in the continuum 
but with a narrow width T/lj p ~ 0.02. The bath temperature is fixed T = 100 and the 
initial temperature of the particle (T ) is varied. We notice that in the case in which the 
temperature of the bath and that of the bare particle are the same, the Boltzmann equation 
predicts no relaxation because the gain and loss processes balance exactly, this is the straight 
line in the graph for bare particle temperature To = 100. The exact solution as well as the 
kinetic and Markovian approximation predict relaxation, the kinetic and the Markovian 
approximations are very close to the exact expression. Analytically we know that the exact, 
Markovian and kinetic will asymptotically approach the Boltzmann result (with very small 
corrections) in this very narrow width case. Obviously the time scales for relaxation and 
the early time dynamics are features not reproduced by the Boltzmann equation and clearly 
a result of off shell effects, since all of the energy conserving detailed balance processes are 
contemplated by the Boltzmann equation. 

Figure [| compares two situations: the left figures correspond to the case of a pole just 
slightly below (dressed particle) and the right figures just slightly above threshold (quasi- 
particle). This case provides for strong renormalization effects because the wave function 
renormalization departs significantly from one. The left figure for g(t) depicts clearly the 
dressing time of the particle, with g(0) = 1 we see that after a short time the asymptotic 
value g{t) « Z cos(u) p t) is achieved. This figure thus reveals two time scales, one associated 
with the oscillation scale of the dressed particle 1/uj p and the other associated with the 
decay to the asymptotic form, this time scale determines the dressing time of the particle 
and for the case under consideration corresponds to just a few oscillations. This dressing 
time scale clearly depends on the details of the spectral density since it determines the early 
time dynamics after the preparation of the initial state. The right figure for g(t) presents 
three different time scales: initially there is the time scale of formation of the quasiparticle, 
very similar to the left figure, the time scale associated with the quasiparticle pole ~ l/oo p 
and finally the time scale associated with the exponential decay. The formation time scale 
and that of exponential decay can only be resolved in the narrow width approximation, in 
this particular example T/u p « 0.005 and the time scales associated with the quasiparticle 
formation from the initial state and exponential relaxation can be resolved. These are clearly 
displayed in fig.|5| where the logarithm of the maxima of g(t) is plotted versus time. In fig|| 
we show the expectation value of the number operator eq.( |2.2| ) for Q = uj p for the same 
values of the parameters as in fig.§ (left figure corresponds to pole below threshold and right 
figure to the pole above threshold) and equal particle and bath temperature T = T = 10. 
Whereas the Boltzmann equation predicts again no relaxation, in the left figure because the 
damping rate vanishes and in the right figure because the on-shell gain and loss processes 
balance each other, the exact and quantum kinetics description of relaxation both predict 
non trivial evolution of the dressed particle and quasiparticle distribution functions respec- 
tively. The left figure shows that whereas the quantum kinetic and Markovian evolution 
are not too different from the exact, asymptotically all of them depart significantly from 
Boltzmann. The early time dynamics predicted by the Markovian and quantum kinetics are 
very close to the exact expression. In the right figure, corresponding to a narrow resonance 
we see that asymptotically the quantum kinetic and Markovian evolution asymptotically 
approach the Boltzmann result but obviously the early and intermediate time dynamics is 
remarkably different. Furthermore the exact result reaches an asymptotic value that is very 
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different from Boltzmann, as a result of the strong quasiparticle renormalization effects, with 
Z departing significantly from one [see fig. |J . Despite the fact that the resonance is rather 
narrow, its proximity to threshold results in strong off-shell effects. 

Figure [7| is perhaps one of the most illuminating. The parameters are the same as for the 
right part of fig.^, i.e. with the quasiparticle pole in the continuum and close to threshold, 
the bath temperature is T = 10 and the particle is initially at zero temperature. In the left 
figure we plot the Boltzmann, exact, quantum kinetic and Markovian evolutions respectively 
for the expectation value of the number operator ( |2.2|) for Q = lu p , whereas the right figure 
corresponds to dividing by Z the results of the exact, quantum kinetics and Markovian 
evolutions to make contact with the quasiparticle number operator ( |6.8| ). This figure clearly 
shows that the Boltzmann approximation coarse grains over the early time behavior and 
completely misses the formation time scales and the early details of relaxation. 

Finally, fig.|| presents the evolution of the quasiparticle distribution for a case of a strong 
coupling regime resulting in a wide resonance: 77 = 5 ; u t h = 5.0 ; u c = 40 ; u p = 
9.58 ; Z = 0.982 for a bath temperature T = 200 and zero initial particle temperature with 
T/u p ~ 0.1. The insert in the figure displays the early time behavior. We see in this figure 
that whereas the early time behavior is similar for the exact and approximate evolutions, 
which is a consequence of zero initial temperature for the particle, at times of the order of the 
relaxation time there is a dramatic departure. Furthermore the Boltzmann approximation 
predicts a very different early time evolution because it coarse grains over the formation 
time of the quasiparticle. 

Whereas the quantum kinetic evolution and its Markovian approximation track very 
closely the Boltzmann, the exact evolution is approximately 15% smaller resulting in a 
smaller population of quasiparticles asymptotically. The departures in the exact result are 
a consequence of off-shell effects associated with a large width of the quasiparticle. 

IX. CONCLUSIONS AND IMPLICATIONS 

The goal of this article is to study the dynamics of thermalization including off-shell 
effects that are not incorporated in a Boltzmann description of kinetics. In particular the 
focus is to assess the validity of the Boltzmann approximation as well as non-Markovian and 
Markovian quantum kinetic descriptions of relaxation and thermalization in a model that 
allows an exact treatment. 

Although the model treated in this article allows an exact solution and therefore provides 
an arena to test the regime of validity of several approximate descriptions of kinetics and 
compare to an exact result, it is obviously not a full quantum many body theory. Specific 
many theory models used to obtain a microscopic description of thermalization and relax- 
ation will certainly contain details that are not captured by the model investigated here. 
However, from the exhaustive analysis in this article we believe that some of the results 
obtained here are fairly robust and trascend any particular model. These are the following: 

• Boltzmann vs. quantum kinetics: A necessary criterion for the validity of a 
Boltzmann approach is that there is a clear and wide separation of time scales between 
the microscopic time scales and the time scales of relaxation. This is typically the 
situation in which quasiparticles correspond to very narrow resonances in the spectral 
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functions and their lifetime is much longer than the typical microscopic scales. If 
perturbation theory is applicable and the quasiparticle resonance is narrow and its 
position is far away from thresholds, then a Boltzmann description is likely to be 
reliable for time scales longer than the formation time of the quasiparticle. When 
there is competition of time scales or the early stages are experimentally accesible a 
full quantum kinetic equation must be obtained. 

Microscopic time scales: In order to determine the microscopic time scales the first 
step is to determine the position of the resonances or quasiparticles, i.e. the quasi- 
particle pole including the medium effects. The bare particle poles do not determine 
the microscopic time scales. Obviously for weakly interacting theories the position of 
the bare and quasiparticle poles will be very close and the microscopic time scales are 
similar. 

Relaxational time scales: An estimate of the relaxational time scale is determined 
by the width of the resonance, T, a wide separation of time scales that would provide a 
necessary condition for the validity of a Boltzmann approximation would require that 
Yjujp « 1. 

Thresholds: Although a wide separation of time scales is a necessary condition for 
the validity of a Boltzmann approach, it is not sufficient. In particular when the 
position of the resonance is too close to threshold, there will be important corrections 
to the long and short time dynamics arising from the behavior of the spectral density at 
threshold. Threshold effects can lead to strong renormalization of the amplitude of the 
quasiparticle pole (wave function renormalization) that results in sizable distortions of 
the equilibrium distributions as compared to the free particle ones. In particular we 
have seen how thermalization is achieved but with large corrections in the quasiparticle 
distribution functions from the usual Bose-Einstein form. The relevance of threshold 
effects can be quantified by the ratio Y / {u p — oj t h)- If this ratio is << 1 then threshold 
effects will be negligible. When this ratio is 0(1) these effects will be important. 

Formation times, Markovian vs. non-Markovian kinetics: The model that we 
have studied allowed us to explore the concept of the formation time of a quasiparticle. 
This concept is simply unavailable within a Boltzmann approach, since the Boltzmann 
equation coarse grains over the formation time scales. This is clearly revealed in figures 
|5] and |7|. Both the non-Markovian quantum kinetics and its Markovian approximation 
include off-shell effects and capture the early time dynamics associated with the for- 
mation of the quasiparticle. The formation time of a quasiparticle becomes relevant 
if the initial state is very far from equilibrium, since in a non-linear evolution, large 
initial departures can result in large corrections in the asymptotic region. It is also im- 
portant if the early time dynamics is resolved experimentally as seems to be the case in 
femtosecond resolved studies of hot carriers in semiconductors. An important message 
learned in this work is that even in strongly coupled cases a non-Markovian quantum 
kinetic description provides a very good approximation to the correct dynamics for 
most of the relevant time scale. A Markovian approximation that is obtained by ex- 
tracting the distribution functions from inside the non-local time integrals, but without 
taking the interval of time to infinity (completed collision) offers a viable description, 
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which is close to the exact evolution and that of the non-Markovian quantum kinetics 
at weak and intermediate couplings. Its main advantage is computational because this 
approximation provides a local update equation. 
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<TG 1. The complex contour used to evaluate g{t) for the cases in which a) the pole is below 
the threshold and b) the pole is above threshold. 
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FIG. 2. The expectation value of the particle occupation number given by eq. (|2.2| ) for Q = oj p 
for the case in which the pole is below threshold for particle temperature To = 200 and bath 
temperature T = 100. The pole is at lo v = 1.95719 and Z = 0.95621. The numerical parameters 
are r\ = 0.85, uj c = 45, u> t h = 5 and ujr = 2. 
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FIG. 3. The expectation value of the particle occupation number given by eq. (|2.2| ) for Q = io v 
for the case in which the pole is above threshold for particle temperatures To = 100 and 200 
and bath temperature T = 100. The pole is at co p = 9.83397 and Z = 0.99631. The numerical 
parameters are r\ = 0.85, u> c = 45, toth = 5 and lor = 10, resulting in T/u> p « 0.02. 
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t t 

FIG. 4. The functions S(u>) and g(t) for the cases in which the pole is just below (left 
column) and just above (right column) the threshold frequency {oJth = 5). In the left column, 
ujr = 5.52, uj p = 4.98083 and Z = 0.63845 while in the right column ujr = 5.65, lu p = 5.07373 and 
Z = 0.69959. The numerical parameters are r\ = 3.0 and u> c = 55, with T/u> p rj 0.005. 




FIG. 5. The logarithm of the maxima of g(t) vs t for r) = 3.0 and tu c = 55, Uth = 5, 
ujr = 5.52, ui p = 4.98083 and Z = 0.63845, corresponding to the right column of fig.|4[ 
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FIG. 6. The expectation value of the number operator eq. (|2.2| ) for the same values of the 
parameters as in figure^) (left figure corresponds to pole below threshold and right figure to the 
pole above threshold) for the case of equal particle and bath temperature Tq = T = 10. The insert 
shows the early time behavior. 
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FIG. 7. The expectation value of the number operator eq.(2.2) for the same values of the 
parameters as for the right column in figure (|j) (quasiparticle pole above threshold), rj 



55, uj t h = 5, lor = 5.65, u p = 5.07373 and Z 



3.0 and 

0.69959. The temperature of the bath is 



T = 10 and zero initial particle temperature (To = 0), the exact, Markovian and kinetic curves 
have been divided by the wave function renormalization Z in the rightmost figure. The insert in 
the left figure shows the early time behavior. 
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FIG. 8. The expectation value of the number operator eq^^) for r\ = 5.0 and uj c = 40, 
uj t h = 5, u! p = 9.58, Z = 0.982 and bath temperature T = 200, corresponding to T/u> p « 0.1. 
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